/*
 * MathFunctions.cpp
 *
 *      Author: Jerome RD Soine and Christoph A Brand
 *      Institution: Schwarz goup, Institute for Theoretical Physics, Heidelberg University, Germany
 */


#include "../../../../include/base/substrate_models/FEM/MathFunctions.h"

/**
 * @brief Return the value of the Levi-Civita Symbol of for a given set of three indeces.
 *
 * This function implements the Levi-Civita Symbol. It is often used to calculate the determinant of matrices
 * or the product of two vectors. The function returns the entry for a given set of indices e_ijk.
 */
double MathFunctions::levi_civita(unsigned int i, unsigned int j, unsigned int k)
{
	if((i==0 && j==1 && k==2) || (i==2 && j==0 && k==1) || (i==1 && j==2 && k==0)) return (1.0);
	else if((i==1 && j==0 && k==2) || (i==0 && j==2 && k==1) || (i==2 && j==1 && k==0)) return (-1.0);
	else return (0.0);
}

/**
 * @brief Returns the value of the Dirac Symbol (or Kronecker Symbol) for a given pair of indices.
 *
 * This function implements the KroneckerSymbol. It corresponds to a unit matrix and return 1 for
 * a pair of to equal indices and 0 otherwise.
 */
double MathFunctions::dirac(unsigned int i, unsigned int j)
{
	if(i==j) return (1.0);
	else return (0.0);
}

/**
 * @brief Calculates the square of a matrix A^2=A*A.
 *
 * Calculates the square of a matrix.
 */
void MathFunctions::square_matrix(double (&matrix) [3][3], double (&squared_matrix) [3][3])
{
  double tmp = 0;
  for(unsigned int i=0;i<3;i++)
    for(unsigned int j=0;j<3;j++)
    {
      tmp=0;
      for(unsigned int k=0;k<3;k++)
	tmp += matrix[i][k]*matrix[j][k];
      squared_matrix[i][j] = tmp;
    }
}

/**
 * @brief Calculates the transposed  of a 3x3 matrix.
 *
 * Transpose function for a matrix. Note, the function is hard coded for a 3x3 matrix.
 */
void MathFunctions::transpose(double (&matrix) [3][3])
{
   double transposed_matrix [3][3];
   for(unsigned int i=0;i<3;i++)
     for(unsigned int j=0;j<3;j++)
       transposed_matrix[i][j] = matrix[j][i];
   for(unsigned int i=0;i<3;i++)
     for(unsigned int j=0;j<3;j++)
	matrix[i][j] = transposed_matrix[i][j];
}

/**
 * @brief Calculates the trace of a 3x3 matrix
 *
 * Calculates the trace of a matrix. Tr(A) = A_ii. Note, the function is hard coded for a 3x3 matrix
 */
void MathFunctions::trace(double (&matrix) [3][3], double &trace_value)
{
  trace_value = 0;
  for(unsigned int i=0;i<3;i++)
    trace_value += matrix[i][i];
}

/**
 * @brief Calculates the determinant of a 3x3 matrix.
 *
 * Calculates the determinant of a (3x3) matrix A. Note, the function is hard coded
 * for a 3x3 matrix and it uses the rule of sarrus.
 */
void MathFunctions::det(double (&matrix) [3][3] , double &det_value)
{
    det_value = 	matrix[0][0]*matrix[1][1]*matrix[2][2]
		+	matrix[0][1]*matrix[1][2]*matrix[2][0]
		+	matrix[0][2]*matrix[1][0]*matrix[2][1]
		-	matrix[2][0]*matrix[1][1]*matrix[0][2]
		-	matrix[2][1]*matrix[1][2]*matrix[0][0]
		-	matrix[2][2]*matrix[1][0]*matrix[0][1];
}

/**
 * @brief Matrix multipication of two 3x3 matrices.
 *
 * Multiplies two 3x3 matrices m1 and m2. Note, the function is hard coded for 3x3 matrices.
 */
void MathFunctions::multiply(double (&m1) [3][3], double (& m2) [3][3], double (&result) [3][3])
{
    double tmp;
    for(unsigned int i=0;i<3;i++)
      for(unsigned int j=0;j<3;j++)
      {
	tmp=0;
	for(unsigned int k=0;k<3;k++)
	  tmp+=m1[i][k]*m2[k][j];
	result[i][j]=tmp;
      }
}


/**
 * @brief Calculates the tensor product of to second order tensors
 *
 * In particular it calculates m1_ij*m2_kl=m3_ijkl.
 */

void MathFunctions::tensor_product(double (&m1) [3][3], double (&m2) [3][3], double (&result) [3][3][3][3])
{
    for(unsigned int i=0;i<3;i++)
      for(unsigned int j=0;j<3;j++)
	for(unsigned int k=0;k<3;k++)
	  for(unsigned int l=0;l<3;l++)
	    result[i][j][k][l] = m1[i][j]*m2[k][l];
}

/**
 * @brief Calculates the deviator of a second rank tensor.
 *
 * The definition of the deviator is m1 = m1 - 1/3*tr(m1)*I and in this way it is implemented.
 */

void MathFunctions::deviator(double (&m1) [3][3], double (&result) [3][3])
{
    double tr_m1;
    double c = 1.0/3;
    trace(m1,tr_m1);

    for(unsigned int i=0;i<3;i++)
      for(unsigned int j=0;j<3;j++)
	result[i][j] = m1[i][j]-c*tr_m1*dirac(i,j);
}


